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In this paper we show how the covariant gauge invariant equations for the evolution of scalar, 
vector and tensor perturbations for a generic /(7?)-gravity theory can be recast in order to exploit 
the power of dynamical system methodology. In this way, recent results describing the dynamics 
of the background FRW model can be easily combined with these equations to reveal important 
details pertaining to the evolution of cosmological models in fourth order gravity. 
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I. INTRODUCTION 

The observational evidence for the accelerated expansion rate of the universe, and the introduction of the concept 
of Dark Energy has put theoretical cosmology into crisis. This is due to the fact that despite an increasing amount 
and quality of data, no model has been proposed thus far that is able to give a completely satisfactory theoretical 
explanation of these observations. 

Among the many different ways to achieve cosmic acceleration, the modification of gravity and in particular higher 
order gravity has recently gained much attention [1-5]. The reason for this popularity is due to the fact that these 
models provide a somewhat more natural explanation of the cosmic acceleration: this effect is due to corrections 
to Einstein gravity which are directly related to the characteristic properties of the gravitational interaction. Most 
investigations of higher order gravity have focused on Fourth Order Gravity (FOG), i.e., on gravitational Lagrangian's 
in which the corrections are at most of order four in the metric and in what follows we will also focus on these models. 

Because the field equations resulting from FOG are extremely complicated, difficult conceptual and technical issues 
arise which need to be resolved in order to uncover the detailed physics of these models. Consequently it is important 
to develop new methods which are able to assist in resolving these problems. Two such approaches, the dynamical 
systems approach to cosmology and the covariant approach to cosmological perturbations, have proved particularly 
useful in this respect. 

The dynamical system approach, first developed by Collins [6] and extensively reviewed in the book by Ellis and 
Wainwright [7] has proved to be an important tool in the understanding of cosmology of /(i?)-gravity models. In fact, 
studying cosmological models using these techniques has the advantage of providing a relatively simple method for 
obtaining exact solutions, which appear as fixed points of the system and obtaining a global picture of the dynamics of 
these models. Consequently, such an analysis allows for an efficient preliminary investigation of these theories, allowing 
one to identify specific models which merit further investigation. In a series of recent papers, a wide range of features 
of f{R) cosmology have been presented, ranging from an analysis of the standard Friedmann-Lemaitre-Robcrtson- 
Walker (FLRW) models [8-10] to a discussion of the properties of the Einstein universe and the isotropization of 
Bianchi models [11-13]. 

The 1+3 covariant approach was developed originally to study the evolution of linear perturbations of FRW models 
in General Relativity in[14-21] with great success. For our purposes this approach has two major advantages. Firstly, 
a specific recasting of the field equations allows one to easily make extensions to a wide range of modified gravity 
theories including Braneworlds [22] and FOG [23, 24]. Secondly the structure of the formalism is such that, unlike 
other approaches , it is possible to keep track of the physical meaning of the equations at any stage of the calculations, 
which can be crucial when one modifies the theory of gravity. 

In particular, in a number of recent papers [23-25] we derived the evolution equations for scalar and tensor pertur- 
bations of a subclass of fourth order theories of gravity characterized by an action which is a general analytic function 
of the Ricci scalar. The results obtained in [23-25] clearly demonstrated that the evolution of scalar perturbations is 
determined by a fourth order differential equation rather than a second order one. This implies that the evolution of 
the density fluctuations contains, in general, four modes rather than two and can give rise to a more complex evolution 
than what is obtained in General Relativity. It was also found that the perturbations depend on the scale for any 
value of the equation of state parameter of standard matter (while in GR the evolution of the dust perturbations are 
not scale dependent) and that there is a characteristic scale-dependent signature in the matter power spectrum [25]. 
This means that in FOG the evolution of super-horizon and sub-horizon perturbations are different. It also turns out 
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that the growth of large density fluctuations can occur also in backgrounds in which the expansion rate is increasing 
in time. This surprising result is strikingly different with what one finds in GR and could provide a strong constraint 
on some FOG theories using the ISW effect. In addition to that the structure of the general fourth order perturbation 
equations and the analysis of scalar perturbations lead to the discovery of a characteristic signature of fourth order 
gravity in the matter power spectrum [24] , the details of which have not seen before in other works in this area. This 
could provide a crucial test for fourth order gravity on cosmological scales. 

The aim of this paper is to combine dynamical systems methods with linear perturbation theory in such a manner 
that one is able to apply directly the results coming from the former to the latter and is able to gain a semi-qualitative 
idea of both the behavior of the background and that of the first order perturbations in a general /(i?)-gravity theory. 
In order to achieve this we will express the coefficients of the perturbation equations in terms of the dynamical system 
variables in such a way that at any fixed point it will correspond a set of perturbation equations and as a consequence 
an evolution law for the linear fluctuations. 

The paper is divided as follow. In Section 2 we give some basic equations. In Section 3 we present a formalism that 
allows one to investigate a large class of /(i?)-gravity models. In Section 4 we present the equations for the evolution 
of the scalar, vector and tensor perturbations in a generic /(i?)-gravity theory written in terms of the dynamical 
system variables. In Section 5 we consider two simple examples. Finally, Section 6 is devoted to the conclusions. 

Unless otherwise specified, natural units {h = c = ks = 9,itG = 1) will be used throughout this paper, Latin indices 
run from to 3. The symbol V represents the usual covariant derivative and d corresponds to partial differentiation. 
We use the — , +, +, + signature and the Riemann tensor is defined by 

R'bcd = W"m,c - Wbcd + W%dW\e - WhcWdf , (1) 

where the W^m are the Christoffel symbols (i.e. symmetric in the lower indices), defined by 

= \9°'^ {9be,d + 9ed,b - 9bd,e) ■ (2) 

The Ricci tensor is obtained by contracting the first and the third indices 

Rab = g'"^Racbd ■ (3) 

Finally the Hilbert-Einstein action in the presence of matter is given by 

(4) 



/ 



II. BASIC EQUATIONS 



In four dimensional homogeneous and isotropic spacetimes i.e., Friedmann Lcmaitre Robertson Walker (FLRW) 
universes, the most general action for fourth order gravity can be written as an analytic function of the Ricci scalar 
only: 



A = J d^x^j[f{R)+Cm] , 



(5) 



where represents the matter contribution. Varying the action with respect to the metric gives the generalization 
of the Einstein equations: 



where / = /(i?), /' = and T*^ 



/'G„b = /' (^Rab - \ 9abR^ = T^b + \9ab (/ - Rf) + ^b^af ' 9abVcV'f' , 

2 5(v^£™) 



(6) 



represents the stress energy tensor of standard matter. 



dR ' ^ Sgab 

These equations reduce to the standard Einstein field equations when f{R) — R. It is crucial for our purposes to be 
able to write (6) in the form 



where = and 



g 

/' 



rpR __ J_ 

ab ^, 



rnvn I rpH rptot 

ab — -I- ab 'T J- ab — ab ! 



^9ab if - Rf) + VbVaf - gabVcV'f' 



(7) 



(8) 
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represent two effective "fluids": the curvature "fluid" (associated witli T^) and the effective matter "fluid" (associated 
with T^) [1, 23]. This step is important because it allows us to treat fourth order gravity as standard Einstein gravity 
plus two "effective" fluids. The details of the conservation properties of these effective fluids have been given in [23] . 
In particular, it has been shown that, no matter how the effective fluids behave, standard matter still follows the 
usual conservation equations = 0. 



III. DYNAMICAL SYSTEM APPROACH AND F(R) GRAVITY 

The DS approach [7] has been used with great success in the analysis of the background evolution of cosmological 
models with fourth order gravity. In [8, 9] it was shown that using the dimensionless variables: 

^3f 3/ 3Mm 

^ /'e' ^ 262' ^ 2/'e2' /'e2' s-^e^ ' ^ ' 

together with the characteristic variable 

Q^bril^ =^ (10) 



d\o%R) Rf" 

and logarithmic time = ] InS'], the cosmological equations are equivalent to the autonomous system: 

dy 
dN 
dz 
dN 
dCt 
dN 
dx 
N 

1 = y + Vl — X — z — X- 



= ey[2x-2y + €li-x + y- z + n-l)+4], 

= sz{Sx-By + z-n + 5) + qey{-x + y-z + n-l), (11) 
= -en{3w-3x + Sy-z + n-2), 

aJS 

dN " 2£x(x-y + i), 



where s = [33]. This system allows one to analyze many interesting fourth order gravity cosmological models 

and leads to the result that some of them present cosmic histories which posses a transient Priedmann phase and 
evolve naturally towards an accelerated (DE-like) expansion phase. 

A detailed analysis of the properties and caveats of this method can be found in [8, 9], here it is important only to 
remind the reader that the solutions associated with the fixed points can be found by substituting the coordinates of 
the fixed points into the system 

6 = 762, j= -'^ + y'-^\ (12) 

. _ (1 + w) 
7 t 

where the subscript "i" stands for the value of a generic quantity at a fixed point. This means that for 7^0, the 
general solutions can be written as 

a = ao(t-to)'/''' , (14) 
Mm = ao{t-to) . (15) 

The expression above gives the solution for the scale factor and the evolution of the energy density for every fixed 
point in which 7 7^ 0. When 7 = the (12) reduces to = which correspond to either a static or a de Sitter 
solution. 

It is important to stress at this point that sometimes the solutions associated with the fixed points are "non physical" 
i.e., they do not satisfy the cosmological equations. Then one might ask how is it possible that, although derived 
from the cosmological equations themselves, our dynamical system equations give non physical solutions. The reason 
for this apparent contradiction needs to be looked for in the very structure of the dynamical system approach. As 
we have seen, the condition to obtain the fixed points of the system (11) -as well as every dynamical system- is to 
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set the first derivative of the dynamical variables to zero (i.e., to set the left hand side of (11) to zero) and solve the 
system obtained. Such a step, in the standard dynamical system theory is usually trivial, however, in our specific 
formulation this step becomes much more subtle. In fact, the requirement of the existence of a fixed point is also 
imposes the requirement that all the variables acquire a constant value (or equivalently that their first derivative with 
respect to the time coordinate is zero). This is equivalent to an additional system of equations that is not necessarily 
satisfied by the solutions of the system. For example, in 

x' = F(x) , (16) 

the condition to obtain the fixed points would be 

F(x) = but also x = const (x' = 0) , (17) 

as mentioned earlier, the second system is usually trivially satisfied. In the formalism above, however, these variables 
are a combination of quantities appearing in the cosmological equations. This means that x = const becomes a 
scit of conditions to be satisfied by all the physical fixed points of the system. In GR, because of the structure of 
the variables [7, 26], these constraints are automatically satisfied. In fact, in this case one has ^gr = 3/Lt™/0^ and 
XGR = 9K/S'^Q'^ which means 

QciJ = const. a , (18) 

XGR = const. 5^ oc KQ'^ , (19) 

so the physical points can either have K = and /x"* oc 0^ or K ^ 0, 5 oc 0. Both the fixed point that one obtains 
(corresponding to Milne and Priedmann solutions) satisfy these criteria. 

This then raises the following question: how do we consider the non physical fixed points? The answer depends 
on the stability. If the fixed point is unstable then, although the solution associated with the fixed point does not 
satisfy the cosmological equations, it can be used to approximate the behavior in neighborhood of the fixed point. 
Orbits nearby this point will bounce on or run away from it, but they will never reach it. Instead in the case in 
which the point is stable, the set of orbits approaching the point will not correspond to any physical evolution for 
the system and the dynamical system approach fails to give an appropriate description of the cosmological evolution. 
This imposes a intrinsic limitation in predictive power of this approach and has to be taken into account to avoid 
incorrect conclusions. 

Another consequence of this limitation is that the structure of the variables characterize the type of solution 
associated with the fixed points and, consequently, the fixed points one derives with a dynamical system formalism 
are not necessarily the complete set of the elementary solutions of the system. This is also important because it 
implies that the absence of a specific cosmic history in a dynamical systems formalism does not necessarily indicate 
that this cosmic history cannot be realized, but only that the specific formalism used is not able to show its presence. 



dN 

dXGR 

dN 



IV. COVARIANT APPROACH TO PERTURBATION THEORY 



The form (7) of the field equations allows us to use directly the covariant gauge invariant approach [14-21] in the 
same way as presented in [23, 24]. Following these references we will choose a frame m™ comoving with standard 
matter {the matter energy frame) [34]. We will also assume that in u™, the standard matter is a barotropic perfect 
fiuid with an equation of state Pm = wfj,m- 

Once the frame has been chosen, the derivation of the kinematical quantities can be obtained in the standard 
manner [14]. In particular the derivative along the matter fiuid flow lines is defined hy X = Ua'V'^X and the projected 
covariant derivative operator orthogonal to is given by Va = h^a^b- With these definitions we can define the key 
kinematic quantities of the cosmological model: the expansion 0, the shear aab, the vorticity coab and the acceleration 
a„ . The general propagation equations for these kinematic variables in any spacetime correspond to the so called 1+3 
covariant equations [14] and are given for completeness in Appendix A. 

The definition of a frame u" also allows us to obtain an irreducible decomposition of the stress energy momentum 
tensor. In a general frame and for a general tensor Tab one obtains: 

Tab = tJ-UaUb + phab + 'iq{aUb) + T^ab , (20) 

where /x and p are the energy density and isotropic pressure, qa is the energy flux and -Wab is the anisotropic pressure. 
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In this way, relative to u^, T*^* can be decomposed as 



tot _ _rptoti^b c _ xm , ti 
% — -'-be "a^ — 9a + 9a 



^tot nntot i-d ^ m i ^ -R 

"^05 — -'cd '^<a'^b> — ^ab '^^abi 



with 



/' 



; — a'" = 



9^ 
/' 



'ab 



"ab 
f 



(21) 
(22) 

(23) 



Since we assume that in our frame standard matter is a perfect fluid, and tt^ are zero, so that the last two 
quantities above also vanish. 

The effective thermodynamical quantities for the curvature "fluid" are 



M = T7 



P = -FJ 



lif- Rf) + f'R + r'R' + l^fR - lf'V'R+ 



'^ab = 



--f"'V''RVaR + f"abV''R 



1 

1 

? 



f'RVaR + f'^aR - ^Qf'^aR 

f"^{a^b)R + f"'^iaRyb)R- <JabR 



(24) 

(25) 
(26) 
(27) 



The twice contracted Bianchi Identities lead to evolution equations for /i^, qj^ and are given in Appendix A. 

Using the quantities defined above, and the equations given in Appendix A, we are able to write both the evolution 
equations for the background and ones of the first order perturbations. As in [23, 24] we will consider a background 
that is homogeneous and isotropic, i.e., a FLRW model. In this background the cosmological equations for a generic 
f{R) read: 



= 3^™ + 3m^ - 



3i? 



+ i(/i" + 3p") + i(Ai^ + 3p^) = , 



fl^ + Q + p™) = 



(28) 

(29) 
(30) 



where R = GK/S"^ is the 3-Ricci scalar, K = 0,±1 and S is the scale factor. 



V. PERTURBATION EQUATIONS AND DSA 

In this section we will rewrite the evolution equations for scalar, vector and tensor perturbations in terms of the 
dynamical system variables. As we will see, the requirement to obtain a closed form for the coefficients will require a 
redefinition of some of the variables in the equations. 

Once this has been done the behavior of the perturbations at a fixed point of the system (11) can be inferred in 
all generality. This can be extremely useful in determining the set of cosmic histories which are compatible with 
observations for a given f{R) model. 



A. Scalar perturbations 

The evolution equations for the scalar perturbations for a generic f{R) theory of gravity in the covariant approach 
were derived in [23, 27] and successively analyzed in detail in [24]. These types of perturbations are characterized by 
the quantities 

S'2 . - - - 

A'" = — V^/x'", Z = S'^V'^e, C = S'^V'^R, n = S^V^R, ^ = S'^VR, (31) 
A* 
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and their evolution equations arc given in Appendix B. In order to express these equations in terms of the dynamical 
system variables we first have to convert them into equations in which the independent variable is N. In addition to 
that one can only obtain a closed form for the coefficients if the curvature variable is redefined as 

R=S^W^ln[f'{R)]. (32) 

Substituting for the new curvature variable and using the definitions (9), the perturbation equations developed in 
harmonics can be written as 

where 

Ai = e(l -3w + z-n), (35) 

Bi = -3(2wz - 3wJC + (1 - w)n) , (36) 

Ci = -3(u> + l)(z-Qy-3/C-17), (37) 

Di = 2,s{w + 1) , (38) 

= -£(3x - 3t/ + 2^ - 20 + 1) , (39) 

Ti=Az- 2Qy - 9/C + (3w - 1)0 , (40) 

w{4.y-^z)-{iw^-Aw + l)n 
^1 = ' (41) 

m = -^-^^^{x-y + z-n + i), (42) 

^ = 2£/C(x-y + l). (43) 

Note that in this form the above equations are such that two different forms of the Lagrangian have the same evolution 

of the scalar perturbations if they both have a fixed point with the same coordinates and and in this fixed point Q is 
the same. As we will see this will allow us to attach a fixed point to a certain evolution law of the perturbations. 



and the form of 1C{N) is given by 



B. Vector Perturbations 



In order to analyze the evolution of the vector perturbations one has to extract the vector parts of the variables in 
the 1+3 equations (A5-A18) and the propagation equations for D^^ Za,TZa and in Appendix B. In our specific case 
some important facts have to be noted. First, looking at (24) the heat flux qa and the anisotropic pressure TVab are 
not independent, i.e., they can be written as functions of the other variables, specifically aab, T^a and Kq. Secondly, 
the variables {Da, ZaTZa, ^a) are one index objects, which means that their purely vector part is obtained by taking 
their curl multiplied by the scale factor. In addition, looking at (31) above one notices that these variables are in fact 
gradients of scalars and we know that curlVo-'i^ = ^UaX. This means that at first order one can write: 

{Daf = -2S^Q{l + w)oJa, (44) 

(Zaf = 2S^eL0a, (45) 

{TZaf = 2S^Rwa = 452 (^lee - + a;„ , (46) 

(5i.)^ = 2S'Ru;a = 45^ (^^|e^ _ ^ X ^ ^ 8^, ^ 8^^ ^ 26(3)) , (47) 

i.e., one can express all these quantities in terms of the vorticity vector. Also one can use the constraints(All- A14) 
to express the remaining quantities aab, Eab, Hab in terms of the vorticity vector. Thus, in the end the only relevant 
equations for the dynamics of the vector perturbations is the vorticity equation (A6): 

e'^ + {2-2,w)iOa = Q, (48) 
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which does not depend on I. This means that the vorticity evolution is independent of the scale. Specifically one has 

J = Jo exp [-eN (2 - 3w)] = JqS-^'^-^'^^ , (49) 

which implies that the vorticity is always decreasing regardless of the features of the action. This is the same result 
that one obtains in GR [17] and shows that f{R) gravity does not affect the evolution of this quantity. However, the 
quantities (44-47) will change behavior according to changes in the background and the form of the action. 



C. Tensor Perturbations 



As already noticed in [28] and shown in Appendix C the only independent equation in the evolution of the tensor 
perturbation is the shear. This happens because the second order equation for aat is closed and from the constraint 
(A12) one obtains Hat = {curla)ab where (curlX)"'' = e'^''^" Vc^^'^d- In addition to that the first order equation for 
cTob (Eq. (CI)) can be used to derive Eab^ because [28] in f{R) gravity (as in the scalar tensor case [32]) the tensor 
component of the anisotropic pressure -Kab can be proven to be proportional to the shear. Thus we are left only with 
the equation for Gab- This equation developed in harmonics and written in terms of the dynamical systems variables 
reads 



^2 =e(2x-2y + ^-0-2), (51) 
B2 = - 2(y - ^ + O - l)x + (y - ^ + 0)2 - 6y + 50 - 9/C + (3w - 2)f2 . (52) 

(53) 

For (50) the same remark given for scalars holds: since the coefficients depend only on the coordinate of the fixed 
points two theories with he same fixed point will have, in the fixed point, the same evolution law for the tensor 
perturbations, like in the case of the scalars. However since the coefficients (51) and (52) do not contain Q such 
occurrences is even more common. 



VI. EXAMPLES 



In this section we will apply the equations defined above to some specific forms of f{R) to illustrate the utility of 
the above approach. 



A. J?"-gravity 



If f{R) = aR" we have a model (often called ii"-gravity) characterized by the action 

L = ^[aR^ + CM], (54) 

which constitutes the simplest possible example of fourth order gravity. In i?"-gravity the characteristic function Q 
is always constant. In particular, we have 

Q = , (55) 

n — 1 

which implies that the variables z and y are not independent, i.e., the phase space of i?"-gravity is contained in the 
subspace y = nz oi the general phase space described by (11). This can be easily seen if one substitutes y = nz into 
(11). Then the equations for y and z turn out to be exactly the same and (11) reduces to : 

dy _ f y~2niy + 2) + 5 , (3 - 2n)x _0 \ 

en (-3w + ( 3 + ) y + 3x - n + 2 ) , (57) 



dN V V 

— = 2£x X + 1 - "i 

dN \ 1 — n 



^^2.,('x + l-T^|, (58) 
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with the constraint 

l + x + y + x-^ = , (59) 

which is equivalent to the one given in [9]. The fixed point with their stabihty and the associated solutions are given 

in Tables I. II and III. In Tabic IV, the long wavelength modes of the solutions for the matter scalar perturbations 
and the tensor perturbations are given. As expected, for the background i2n/3(i+iii) corresponding to the point G, the 
results are the same as the ones already found in [23, 28]. For £ ^ 0, however, one has to use numerical methods to 
obtain the solution of the equations. 

It is interesting to observe the behavior of the matter perturbation modes of the point Here the matter 
perturbations posess a constant mode and the other modes can be growing or decaying depending on the value of the 
parameter n. 

As shown in [10] for some specific values of the parameter n (1.37 ^ n ^ 2) this model has a set of cosmic histories 

characterized by the presence of a transient, decelerated power law expansion that evolves towards an accelerated 
expansion. Using (33) and (50) we are now able to see directly the evolution of the scalar and tensor perturbations 
along these orbits. 

In particular, as the Universe approaches the point !F, for 1.37 ^ n ^ 2 and w = 0, the large scale scalar 
perturbations, which nearby Q have a growing mode (see Figure 1), start dissipating, which is consistent with what 
one would expect in a late time acceleration scenario. The large-scale tensor perturbations instead do not change 
their behavior and keep being dissipated, but at a much faster rate. 

In order to analyze the behavior of the perturbations for smaller scales { £ 0) one needs to integrate the equations 
numerically. This can be done in a relative easy manner and an example of the results obtained in the case of dust 
and £ = 100 are shown in Figures 3 and 4. It is clear that in the point the matter scalar perturbations approach 
a constant value which depends on the initial conditions, while in the point Q the perturbations first have a phase 
of growth and then start to decay which is consistent with what was found in [24]. The same can be done with the 
tensor perturbations, but we find that, as for the £ = 0, case these types of perturbations are dissipated on small 
scales. 




FIG. 1: Plot against n of the real part of the exponents of the long wavelength modes for ii"-gravity in the point Q and in the 
dust case (w = 0) . 



In this case the action reads 



VII. 

L 



THE CASE f{R) = R + aR 



(60) 
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FIG. 2: Plots against n of the real part of the exponents of long wavelength matter perturbation modes for i?"-gravity in the 
point in the case of dust (w = 0). 
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FIG. 3: Plots of the solutions of the equations (33) in the fixed point in the case n w 1.37, w = and I = 100. 
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FIG. 4: Plots of the solutions of the equations (33) in the fixed point Q in the case n w 1.37, w = and ^ = 100 
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TABLE I: Coordinates of the fixed points for the model f{R) = aR". The superscript "*" represents a double solution. The 
point jB is a double solution for n = 0, 2. 

Point Coordinates {x,y,z,0,) X 

A* (0,0,0,0) -1 

B (-1,0,0,0) 

C (-l-3«;,0,0,-l-3w) -1 

V (l-3w,0,0,2-3w) 

£ (2(l-n),2n(n-l),2(n-l),0) 2n(n - 1) - 1 

T- ( 2(n-2) (4n-S)n 4w-S n\ n 

•' \ 2n-l ' 2n^-3ri+l' 2n^-3n+l'"^ " 



Q 


i- 




-1) 4n-:i(ir-l) 
























TABLE II: Solutions associated to the fixed points of the model f{R) = aR". 



Point 




Scale Factor 




A 






Mm = 


B 


S = 


So{t - tof^ (only for n = 3/2) 


Mm = 


C 




S = So{t - to) 


Mm = 


V 


s = 


Soit - to)^^^ (only for n = 3/2) 


Mm = 


£ 




/^=2;i^ if 

\ S = Sot if K = 


Mm = 






(l-n)(2n-l) 
S = Sot n-2 


Mm = 


G 




S' = S'of3(i+™) 


Mm — Mm,0^ 



l^^^o = (-l)"3-"22"-in"(l + w)-^"x 
(4n - 3(1 + w))"-^[2n^(4 + 3w) - n(13 + 9w) + 3(1 + w)] 



This theory has gained much popularity as a fourth order gravity model within the context of both inflation and dark 
energy [3, 29-31]. The characteristic function Q is : 



Q = 



n{z - y) ' 

and substituting this relation into the system of equations (11) we obtain 



dy 
dN 
dz 

m 

dQ 
dN 
dx 
dN 



vi—x + y — z + Q — 1) 
2x-2y+ — ^ ^ + 4 

ny — nz 

(-X + y - ^ + - 1)2/2 



with the constraint 



ny — nz 

-£0(3w - 3x + 3y - ^ + f2 - 2) , 

2ex(x - y + 1) , 



1 + a; - y + ^; + X = , 



+ ^(3x - 3y + z - n + 5) 



(61) 

(62) 
(63) 
(64) 

(65) 

(66) 



The fixed points and their stability for the phase space of (62) are shown in Tables V VI and VII. In Table VIII we 

show the exponent of the long wavelength scalar and tensor perturbation modes. 

The first point to note is that there exists a fixed point that corresponds to a transient Friedmann type behavior just 
as in the f{R) = -R" case for exactly the same values of the dynamical system variables. It is therefore not surprising 
that the long wavelength scalar perturbations modes for these fixed points have the same solutions in both theories. 
At first glance this may seem to contradict the results in [24] , where it was found that the scalar perturbations at the 
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TABLE III: Stability of the fixed points for i?"-gravity with matter. We consider here only dust, radiation. The term "spiral"*"" 
has been used for pure attractive focus- nodes and the term "saddle- focus" for unstable focus- nodes. 



n < 1(1 - V3) i(l - V3) < n < < n < 1/2 1/2 < n < 1 



A 


saddle 


saddle 


saddle 


saddle 


B 


repellor 


repellor 


repellor 


repellor 


C 


saddle 


saddle 


saddle 


saddle 


V 


saddle 


saddle 


saddle 


saddle 


£ 


saddle 


attractor 


spiral 


spiral 




attractor 


saddle 


saddle 


attractor 



1< n < 5/4 5/4 < n < 4/3 4/3 < n < ^(1 + V3) n>\{l + ^3) 



A 


saddle 


saddle 


saddle 


saddle 


B 


saddle 


repellor 


repellor 


repellor 


C 


saddle 


saddle 


saddle 


saddle 


V 


saddle 


saddle 


saddle 


saddle 


S 


spiral 


spiral 


attractor 


saddle 


JF 


repellor 


saddle 


saddle 


attractor 



g n< 0.33 0.33 < n < 0.35 0.35 < n < 0.37 0.37 < n < 0.71 0.71 < n < 1 
w = saddle saddle-focus saddle-focus saddle-focus saddle 
w = 1/3 saddle saddle saddle- focus saddle- focus saddle- focus 

1 < n < 1.220 1.220 < n < 1.223 1.223 < n < 1.224 1.224 < n < 1.28 

w = Q saddle-focus saddle-focus saddle-focus saddle-focus 

w = 1/3 saddle- focus saddle- focus saddle- focus saddle- focus 

1.28 < w < 1.32 1.32 < w < 1.47 1.47 <n< 1.50 n > 1.50 

w = saddle-focus saddle saddle saddle 

w = 1/3 saddle saddle saddle saddle 



fixed point C depend on the value of a. However, this discrepancy can be resolved when one considers the structure 
of our dynamical system formalism and in particular the conditions (17). These additional equations impose further 
constraints on our system. This can be seen more clearly if we consider the fixed point £ in the case of dust {w = 0) 
in more detail. Using the definition of the dynamical systems variable, x and our choice of f{R) we find that the 
following constraint must be obeyed at C 

3(n - 1) S an'^RR^''-'^'^ 
x = — 5^ ^ — = -1. (67) 

If we then use the fact that 




^=6 o + + . 5(i) = 5o(i-to)^ (68) 



and that we require that the constraint above is satisfied for all time, it is trivial to show that (67) is only satisfied 
if one requires a — > oo. This means that supposing that the point £ is associated with a solution of the cosmological 
equations, means that the function f{R) we are dealing with is very close to RT" . As a consequence, when we insert 
constant values of the dynamical system variables in the coefficient of the perturbation equations wc arc in fact 
imposing that we are dealing with a theory which is essentially i?"-gravity. Thus, naturally, we recover exactly the 
same results of the previous section. This can be also seen from the fact that if we substitute the coordinates of the 
fixed point in the definition of Q we obtain 



(69) 
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TABLE IV: Exponent of the modes for scalar and tensor in the fixed points of the system (56) . 



Point Scalar Modes Exponents 



Tensor Modes Exponents 



A 
B 
C 

V 
£ 



{-2e,0} 

{0,£,e(^-l)} 

{-2e, 0, e(3w - 1), e(3M; + 1)} 



{-£(2 + V3),-£(2- V3)} 

{-£,0} 

{§ (3w - V9w2 - 6w + 9 - 3) , 
f (3w + V9w2 _ 6u; + 9 - 3) } 

{0,.,e(i^-l))} {0,£(%^-2)} 

|o, £(3m; - 1), £ [n - 2 \/3\/n(3n - 4)1 , 

, e [n - 2 + V3 V'^(3n - 4)j j {e (n - 3 - V6 - in^) , e (n - 3 + V6 - 3n2) } 

{0,e(8n + 2+;^),£[4(n + l) + ;A_], 

{£n(4+^),8e(n + l + ;j^)} 



_g I 3(n(2n-3Hl)"' + l| | 

4n \ / 'Invj 
'Mw + l) ) '£ \ m + l 
3(iu + l) + 3n((2n-3)m-l) + \/A' 

6(n-l)(tu + l) 
3(TO+l) + 3n((2n-3)u)-l)-VA 
6{n-l){u) + l) 



/l 4n \ 2(n-l)m-2\ 

\^ 3(m + l) j -m+l J 



A = (n - l)(4(3w + 9,yrf - 4(3w(18w + 55) + 152)n^ + 3(w + l)(87w + 139)n - 81(w + 1)^) 



which is exactly the Q of i?"-gravity. Since Q is the only parameter that differentiate the theory of gravity in the 
coefficients of the perturbation equations (as well as the dynamical system) we clearly expect that after substituting 
the coordinates of the fixed point into the c;ocfficicnts we obtain the same results of i?"-gravity. 

This seems to suggest that in some sense the fixed points carry information about the theory of gravity other than 
the background evolution, so that a fixed point represents a physical solution only of a specific form of the Lagrangian. 
Hence if the phase space of a generic theory of gravity posses that fixed point it means that there can exist regimes 
in which this theory can be approximated by a Lagrangian for which that background is a "physical" solution. 

Furthermore, from what was said above, one can also conclude that somehow the evolution of the perturbations 
is attached to the fixed point in such a way that regardless of the theory, the evolution of scalar perturbations is 
determined only by the fixed point. In other words we are somehow obtaining the "fixed point" of the perturbation 
theory which corresponds to the fixed point of the dynamical systems approach. Such a fixed point is an approximation 
of the real behavior of the equations. This means that we expect the scalar perturbations around this fixed point to 
be approximated by the results we found, which is consistent with what we obtained in [24] . 



TABLE V: Coordinate of the finite fixed points for R + ai?" gravity. 



Point Coordinates (a;,j/,0, f2) X 

A (0,0,0,0) -1 

B (-1,0,0,0) 

C (-l-3w,0,0,-l-3w) -1 

V (l-3w,0,0,2-3w) 

£ (0,2,1,0) 

^ (2,0,-2,0) -1 

g (4,0,-5,0) 

H (2(l-n),2n(n- l),2(l-n),0) 2n(n - 1) - 1 

7" ( 2(71-2) (5-4n)n 5 -4k r)^ r. 
\ 2n — 1 ' 2n^ — 3n+ 1 ' 2n^-3 

/ 3(n-l)(ro+l) -4n.+3u)+3 -4n+3TO+3 -2(_3w+4)n'^ +(9w+13)n-3(w+l) \ „ 

\ n ' 2n ' 2^ ' 2^ J ^ 



VIII. CONCLUSIONS 



In this paper we have discussed the connection between the dynamical system approach and the covariant-gauge 
invariant theory of perturbations, presenting a method to calculate directly the evolution of the scalar, vector and 
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TABLE VI: Solutions associated to the fixed points of R + aR^. The solutions are physical only in the intervals of p mentioned 
in the last column. 



Point 



A 
B 
C 
V 

S* 

JF 
G 

n 

T 
L 



Scale Factor 



S = S^{t- to) 
S = So{t- tof" 

S = So{t- to) 
S = So{t- tof^"^ 

S = So, 

S = So exp [±2V3aT(2 - 3n)''{t - to)] , 

^ ~ 2{l-n) 

S=it- to) 
S = So{t- tof^ 
S= -v/l - 2n(n - 1) (t - to) 

S = So{t- to) 
S = So{t-to)T!^ 



Energy Density 



H = HmOt 
















3(211 



Physical 



-3.. + l)(ra + l) 



fJ-m = fJ-m,,o{t — tp)^'' 



n > 1 

n > 1 

n > 1 

n > 1 

n > 
n < |,a > V 
n > |, Q < 

n > 1 

n > 1 
l<n> i + ^ 

n = |,/tim,0 = 

non physical 



TABLE VIL The stability of the fixed points in the model R + aR". The quantities Bi related to the fixed point C, represent 
some non fractional numerical values {Bi w 1.220, Bi w 1.224, B3 w 1.470). 

Point Stability 



B 


saddle 




B 


] 


repellor 


< w < 2/3 




saddle 


otherwise 


C 


saddle 




V 


< 


repellor 


2/3 < w < 1 




saddle 


otherwise 



G 

n 



attractor 
spiral"*" 
saddle 
saddle 
saddle 

attractor i(l 
spiral"*" 
saddle 
attractor 



repeller 

saddle 



§ < n < 2 
0<n<i 
otherwise 



- 73) < n < 
< n < 1 
otherwise 
n< i(l-^)Un>2, 
l<n< |,(w = 0,1/3), 
l<n< ^(ll + ^),(w = l) 
otherwise. 





' TO = 0,1/3 




saddle. 


' 1 


w = 1 


1 repellor 


Si < n < S2 U S3 < n< |, 






1 saddle 


otherwise 



tensor perturbations at a fixed point of tlie piiase space of a generic /(i?)-gravity theory. Within the limitations of 
the dynamical system approach one is then able to obtain an idea of the evolution of the perturbations in any f{R) 
model. 

Because of the non-linearity and the peculiar structure of the dynamical system formalism, the concept of fixed 
points in /(i?)-gravity is more subtle than the one of GR. In particular, one can have fixed points which do not 
correspond to solutions which satisfy the cosmological equations. This is due to the fact that the fixed point conditions 
(17) constitute additional constraints that can be incompatible with the cosmological equations and that the exact 
solution associated to the fixed points seems to depend more on the definition of the variables than the model itself. 
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TABLE VIII: Exponent of the modes for scalar and tensor in the fixed points of the system (62). 



Point Scalar Modes Exponents 



Tensor Modes Exponents 



{-2e,0} 

{-2e, 0, e(3w - 1), e(3M; + 1)} 
{0,e,£(^-l))} 

-3n-^/25n-32y/n -3n+\'25n-32y/n 
2n 2n 

{-2s, -2s, 0, -e(l - 3w)} 

{-5e,0,2e,e(3w)-2)} 
|o, e{3w - 1), e - 2 - V3y^n{Sn - 4)] , 

• 2 + V3 Vn(3n - 4)] | 

I {0,£(8n + 2+;^),£[4(n + l) + ;^], 



A 
B 
C 

V 
£ 

G 
H 



} 



{-£(2 + V3),-£(2- V3)} 

{-£,0} 

{§ (3w - V9w2 - 6w + 9 - 3) , 
f (3w + V9w2 _ 6u; + 9 - 3) } 

{0,£(¥-2)} 
{0,0} 

{e(-3- V6),e(-3 + v/6)} 
{-le,0} 




3(n(2n.-3) + l)m 



n-2 
4n 



+ 1 



]} 



3(it 



3(^u+l)+3n((2n-3)w-l)+^/A ] 
6(n-l)(ii;+l) J 
3(-m+l)+3n((2n-3)w-l)-V^ 1 
6(n-l)(w+l) I 



i 



{e (n - 3 - V6 - 3^2) , e (n - 3 + V6 - 3^2) } 
{e^ (4+Tr^),8£ (n + l + ;^)} 



, 2(w-1)tu-2 
' w+l 



A = {n- l)(4(3w + 8)^71=* - 4(3w(18w + 55) + 152)n^ + 3(w + l)(87w + 139)n - 81(u; + 1)^) 



This has profound implications on the interpretation of the results of the dynaniieal system approach. Specifically 
it suggests that this kind of approach offers only partial information on the actual evolution of the cosmology in this 
framework and that this information depends on the formalism used. In particular (i) one might not be able to see 
all the fixed points in the phase space because of the form of the variables and (ii) some of these fixed points might 
not correspond to actual solutions of the system with obvious problems in the interpretations of the orbits which 
have these points as attractors. As a consequence one has to be very careful in using these tools to derive general 
conclusions about the dynamics for these cosmological models. 

These peculiarities also have consequences on the results of the perturbations equations. In particular, we found 
that the evolution laws obtained using this form of the equations do not in general coincide with those that one 
would obtain by simply subsittuting in the background corresponding to the fixed point in the dynamical system 
equations. This can be explained when one considers that the additional conditions associated with the fixed points, 
further constrain the perturbation equations and therefore lead to different results. On the other hand such a fact also 
implies that one can associate an evolution law for first order perturbations to a specific fixed point, which in some 
sense may be considered as a "fixed point" for the perturbation theory. Thus one can use the results of the equations 
given above in the same way in which one used the solution at the fixed points: gaining qualitative information about 
the behavior of the perturbations along an orbit. Such a feature, confirmed by the direct calculations presented in 
[24], can help with the understanding of the behavior of the perturbations in models for which a direct numerical 
integration is to complex or to resource intensive to be performed. It is also worth noticing that since the peculiarities 
in the predictions of the equations presented above derive ultimately from the additional condition on the fixed points, 
they only apply when one deals with fixed points and will not be present when one considers the dynamical system 
variables as functions of time. 

In conclusion, in spite of these difficulties the dynamical system approach, when combined with the covariant gauge 
invariant formalism, represents an extremely powerful tool for the study of /(i?)-gravity cosmological models. We 
believe that a careful use of these methods could be invaluable in determining the relevance of these models on 
cosmological scales and to constrain them using current and future observations. 
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APPENDIX A: GENERAL PROPAGATION AND CONSTRAINT EQUATIONS OF THE 1+3 

COVARIANT FORMALISM. 

In this appendix wc will write the propagation and constraint equations for the kineniatical and thcrniodynamical 
quantities of the 1+3 formalism for a generic /(i?)-theory of gravity. Chosen the frame Ua the kinematical quantities 
are defined as 

e = V"Ua, fTab = V(aW6), Wo6 = VfaUfc] , Ua = UcVUa , (Al) 

and the thermodynamical ones as in (21) and (22) with the condition 

Qa = q{a) ) ^ab = T^(ab) ■ (A2) 

As usual angle brackets applied to a vector denote the projection of this vector on the tangent 3-spaces 

V(^a) = K'Vb . (A3) 
Instead when applied to a tensor they denote the projected, symmetric and trace free part of this object 

W(a6> = [hia'K/ - ^h-'hab] • (A4) 

Expansion propagation (generalized Raychaudhuri equation) : 

e + + aab'j"'' - 2a;„a;" - V"ua + iiau" + ^(/i'" + 3p™) = -i(/x« + 3p«) . (A5) 
Vorticity propagation: 

W(a> + fQWa + jCUrlUa - CTafcO;^ = . (A6) 

Shear propagation: 

O'(afc) + §0O-ab + Eab " ^ (aUb) + CTc{a<^bf + <^{a<^b) " "(a^fc) = ^T^ab ■ (^7) 

Gravito-electric propagation: 

E(ab) + QE„,b - curlHab + i(A" +P")'7afc " 2w^£,rf(„i/fc) - Sa^^aEb)'' + iv'^e^diaEb)'' 

- +P )^ab - 2^{ab) - 2^{a%) ' 6*^^a6 " 2^ (o.^b)c " 2^ ^c{aT^b)d ■ (A8) 

Gravito-magnetic propagation: 

H(ab) + ^Hab + CUrl_Eah — 5ac{aHb)'^ + '^'^^ cd{aHb)'^ + '^u''Scd{aEb)'^ 

= bmln^b - i^iaQb) + ^'ia^bydd ■ (A9) 

Vorticity constraint: 

V^Wa - il^Ua = . (AlO) 

Shear constraint: 

VVafe - curl - fVaO + 2[u;, u]a = -Q^ • (All) 

Gravito-magnetic constraint: 

CUrl(To6 + V(aa;6) - Hab + 2U(aW6) = . (A12) 

Gravito-electric divergence: 

V^K6 - I V„/i™ - [<7, H]a + m^buj'' = lallb - 1 [w. qX - h^'^ab + i VaM^" - |eg^ . (A13) 
Gravito-magnetic divergence: 

y'Hab - ir+Pn^a + W,E]a - ^Eabuj' = -^curlgf + (/x« +p«)u;„ - \[a,nX " hab^" ■ (AM) 
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Standard Matter Conservation (twice contracted Bianchi identities for standard matter) 

= -e(M'"+p") , (A15) 

V>™ = -(/x'" + p™) u« . (A16) 
Curvature fluid Conservation (twice contracted Bianchi identities for the curvature fluid ) 

+ W^qR = _ e (/x« + p«) - 2 («»gf ) - (a^^Trfa) + /x™^ , (A17) 

Qfa) + VaP^ + V^S, = -iQq^-aa' - {n"" + P^) < - - vl' ^6 + M™^^^^ • (A18) 

In the equations above the spatial curl of a vector and a tensor is 

(curl Xy = e«^^ VbXc , (curl Xy^ = e"'^'^" VcX% , (A19) 

respectively, where Cabc = u'^Vabcd is the spatial volume. Finally uJa = ^Sa'"^^bc and the covariant tensor commutator 
is 

[W,Z]a=SacdW%Z'"'. 

The 1+3 equations above arc completely equivalent to the Einstein equation and govern the dynamics of the matter 
and gravitational fields in fourth order gravity. As we will see the new source terms in their R.H.S. will modify the 
evolution of the perturbations in a non-trivial way. The standard OR equations are obtained by setting f{R) = R 
which corresponds to setting all these sources to zero. 



APPENDIX B: THE COVARIANT GAUGE INVARIANT PERTURBATION EQUATIONS FOR 

SCALARS 

In the following we give the equations for the evolution of the gradient variables as done in [23, 24] 
They read 



Rf" 


29 


f 


" IT 


1 





3(u>- l)(3w + 2)^'" 2we2 + 3w(m^ + 3^-^) 2w K 
6(w + 1) Y ^ 6(w+ 1) ^ w + l'S^ 



2/' " 



f 



7^a 



w 



5Ra 



= 5?a- 



-(e + 2i?^^3f?a-^Z„- 



(3«; - 1) /z™ „ w , E ^?^ /' w 

Q 777 + 3— -t p'^ + m'' 777 + o/ 
3 /" u> + 1 /" 3(w + l) 

-{l+w)^-y-— + —R +0—R--0 +3p ) + Rjjr 

1 f 1 /;'"!• f" 



R{e-3R^^ 



(Bl) 
(B2) 

(B3) 

(B4) 



a 



6/' ' 2^-" ' ^' f 
together with the constraint 

2Rf"^ 



/' 



(B5) 



^ ^ ^ |e + Za - 2^p™ + [^RQyr - 772 (/ - 2m" + 2i?e/" - 4^ j J + - ^y^Ha = o . 

(B6) 



29/'' 2/"--, 
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The propagation equation for the variable Ca is 



Ca = 



6K^ 

5*2 e 



5*2 e 



Ca 



Qjj^ 2(3a;)9 
~e~ ~ 3(a; + 1) 



+ 



6it!7"2 (2 (92 - 3^1^) f - 3/) /" 6i?7(3) 

~l~ fin ~l~ 



e/'2 



/' 



7^a 



4u;52e 252/" 252(e/"-3E/(3)) 

3(0; + 1) " + /' " 3/' 



(B7) 



If we focus on the evolution of scalar part of these variables, which is associated with the spherically symmetric 
collapse terms and can be extracted taking the divergence of the definitions (Bl), these equations become 



wOA^ - (1 + w)Z . 



l_llfl_r[^ 
2 2 /' /' /' /' 

w 



(w - l)(3w + 2) /x™ _^ 2^92 + 3w(/i^ + 3p^) 
2iw + l) 



6(w; + 1) 



3f? + 



w + 1 



SR = - e + 2i? 



f" J 



^-RZ 



f 

[ZW — 1) jJL 



. f(3) 



w + l 



-V2A. 



J Til 



+ 



w 



f" 3(u. + 1) 



R 



(B8) 

(B9) 
(BIO) 

(Bll) 



C = K' 



18/"7^ 18A„ 



s'2e/' 529 



K 



C + A„ 



2(«;-l)9 6^*^ 
w + l ~9~ 



6^9/7(3) - /" (3/ - 2 (92 - 3At^) /' + 6R&f") 
J }^ 1-7? 



iws^e ^ , 2sy"^ 252 (9/" - 3i?/(3)) 



3(w; + l) 



/' 



3/' 



together with the constraint 



29 f" 2f" ~ 
7e+— ^JR- -^V27e = 0. (B12) 



Traditionally the analysis of the perturbation equations is simplified by using a harmonic decomposition. In the 
1+3 formalism this can be done by developing the scalar quantities defined above using the eigenfunctions of the 
Laplace-Beltrami operator [16]: 



(B13) 
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where I = 2-!tS/X is the wavenumber and Q = 0. Developing (31) in terms of Q, (B8-B12) reduce to 

weA(^) - (1 + , 

'"(w- l)(3w + 2)^'" , 2u;e2 + 3w(^« + 3p«)] ^(^) 



m 



Rf" 29 



2iw + l) f 



+ 



1 /" 1 / /" /" A^" 



Re 



f" 



w + l 



/' 

(3'u; — 1) /i 



6(w + 1) 

/(3) 



9/" 



/' 



/' 



+ 



J2 n f J(4) . y(3) j(3) 



52 

18/"7^ 18A„ 



/" 3(u; + l) 
R 
"3 



aw + 



6^1 srW 

e/' 



+ 



529/' s'29 

652i?9/7(3) - 6£2/"/' + 52 J-// |^3j _ 2 (92 - 3/x«) /' + 6i?9/") 



+ g2 



95'2(/')2 

4^^2e ^^^^ ^ 252/-^,, 2^2(er_3^^(3)) 



+ 



3(w + l) 



/' 



3/' 



(B14) 

(B15) 
(B16) 

(B17) 



(B18) 



3^+ /' i /' - 



■ f f " / ■ \ f " ^2 ■ 



/' 



(B19) 



which is a system of ordinary differential equations. This system takes a more manageable form if we reduce it to a 
pair of second order equations: 



= -{w + l) 



72 \^ Rf"' 


A(^)- 




"to 



'52 



-u;(3p« + M«) 



2wRQf" (3^2 - 1) /i" 



a(^) 



(2 f„ . ^„ 



2i?9 



/(3) 



f ' V - ■ " y /'2 /' 

2 /<" 9 //"^ 1 



/' 



/"^W + (9/" + 2E/(3))7tW- i^. 

f f//2 



7^(^) = - 



2/' 

-(3w;- 1)m'"+ 



+ 



(B20) 



(B21) 



APPENDIX C: THE COVARIANT GAUGE INVARIANT PERTURBATION EQUATIONS FOR 

TENSORS 



These equations were derived in [28] here we sketch briefly their derivation. If one focus only on tensor perturbations 
the 1+3 equations reduce to 



2 1 

(^ab + 2 0fo6 + Eab — 'j^^ab = 



Hab + Hab<d + (curl£;)ab - ^(curl7r)afe = , 



(CI) 
(C2) 
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Eab + EabQ- {cml H)ab + ^ {iJ.*"* + P*"*) (Tab + TTab + ^T^ab = , 

together with the conditions 

VhH^'' = , VbE^'' = , Hab = (curl aU ■ 
Taking the time derivative of the above equations wc obtain 

(Tab - VV + ^ e (T„6 + Q + - aab = 7r„6 + ^ 6 TTab , 



(C3) 



(C4) 



(C5) 



iiab - y'Hab + ^ e TJab + I - 3p*°') Hab = (Curl7^)„fc + I e (curl n)ab 



(CG) 



Eab - ^^Eab +i:QEab+- - 3p) Eab " ^6 (/Z*°* + (1 + 3 <7ab 



l^Tab - ^ VV„6 + ^ e + i (e^ - 7r„6 



(C7) 



where the sound speed is defined as = p*"* / and we have used the Raychaudhuri equation (Eq. (A5)), the 
matter conservation equation (Eq. (A15)) and the commutator identity 



(curlX)„b = (curlX)„, + -(curlX) 6 . 
Using trace-free symmetric tensor eigenfunctions of the spatial the Laplace-Beltrami operator defined by: 



(C8) 



^^Qab = oQab 



(C9) 



where £ = 21^8/ \ is the wavenumber and Qab = one obtains 
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rII I ifW = 0, (Cll) 



(C12) 



where we have used (CI) as an equation for Eab and we have used the fact that the tensor part of the anisotropic 
pressure is Wab = jT(TabR- It is clear that the first two equations are closed and the third is only a constraint. 
Therefore the only relevant equation in our case is (CIO). 
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